% This is the main file to run for "The Climate Risk Premium" (Lemoine,
% JAERE, 2020).

% Needs the compecon toolbox from Miranda and Fackler (2002)

clear all;
close all;


%% User options

% main directory
dirmain = 'C:\yourdirectory';
customdir = '2020-01-23';  % header for subdirectory under dirmain in which will save workspace and any plots

% Numerical parameters
Params.montecarlo_draws = 1e3; % default: 1e3; for consumption volatility and, if in that mode, for simulations
Params.quadrature_nodes = 8; % default: 8; quadrature nodes for climate sensitivity and, with limited flexibility, for damages (which are hard-coded)
Params.horizon = 200; % horizon, in years; default: 200
Params.years_for_scc = 1; % number of years in which will calculate scc; default: 1
Params.timestep = 1; % timestep for discretizing stochastic processes, in years, should be a fraction <= 1; default: 1

% Determine whether want to loop over parameters or settings
option_otherloop = 'off'; % 'off', 'forms_of_unc', 'cons_vol', 'temp_vol', 'cs_var', 'damage_var', 'cs_and_damage_var', 'eta', 'eta_with_stern', 'emint', 'damage_exp', 'damage_exp_with_stern', 'cons_growth', 'discount', 'montecarlo_cs', 'montecarlo_damage'

% Types of uncertainty to include
option_uncertain_climsens = 'on'; % 'on' or 'off'; 'on' means have uncertain climate sensitivity
option_uncertain_damage = 'on'; % 'on' or 'off'; 'on' means have uncertain damages
option_consvolatility = 'on'; % 'on' or 'off'; 'on' means have nonzero volatility for consumption
option_tempvolatility = 'off'; % 'on' or 'off'; 'on' means have nonzero volatility for temperature

% Robustness options
Params.perturbation = 0; % =0 (default): use analytic scc; >0: finite-difference calculations, done by adding this many units of CO2 (Gt C) to time 0 and calculating change in welfare
option_damagetype = 'base'; % 'base', 'level', or 'nordhaus'
option_ocean = 0; % =0 (default): one-layer model without distinct ocean sink; =1: include deep ocean for two-layer model
Params.do_ez = 0; % =1: do Esptein-Zin value recursion; =0 (default): don't

% Other options
option_discount = 'nordhaus';  % 'nordhaus', 'stern', 'halfnordhaus', or 'doublenordhaus'
option_nodamages = 'off'; % 'on' or 'off' (default); 'on' means that warming does not affect consumption

option_dohpc = 0; % =1: using high performance computing facility; =0: using personal computer


%% Recursive utility options - do not need to be adjusted

% can adjust preference parameters below: Params.risk and Params.time

Params.learn_ez = 1; % =1 (default): learn climate sensitivity and damage parameter after first instant; =0: never learn them; only relevant if Params.do_ez=1

% Computational parameters for Epstein-Zin
Params.basiscoeffs = [250 8 8 8]; % 4x1 vector defining order of approximation, ordered consumption p.c. then M1 then M2 then T; default: [250 8 8 8]
Params.nodes = Params.basiscoeffs; % 4x1 vector defining order of approximation, ordered consumption p.c. then M1 then M2 then T; default: Params.basiscoeffs
Params.basistype = 'lin'; % 'cheb', 'spli', or 'lin' (default); type of basis functions to use for approximation
Params.nodetype = Params.basistype; % 'cheb', 'spli', or 'lin'; default: Params.basistype; type of collocation nodes
Params.quadrature_nodes_consvol = 5; % default: 5; number of quadrature nodes to use when taking expectations over consumption shocks
Params.quadrature_nodes_consvol_statespace = 8; % default: 8; number of quadrature nodes to use when simulating consumption shocks in order to bound the state space; should be > Params.quadrature_nodes_consvol 
Params.resetjumps = 0; % =1: reset all jumps past grid back to grid; =0 (default): don't
if strcmp(Params.nodetype,'cheb') % cheb works better without the log transform, but state space so big that starting state often not contained in it
    Params.uselogstates = 0;
elseif strcmp(Params.nodetype,'spli')
    Params.uselogstates = 1;
else
    Params.uselogstates = 1;
end
if option_dohpc~=1
    Params.workers = 5; % number of workers to use when parallelizing nodes, should be weakly positive integer
else
    Params.workers = 28; % number of workers to use when parallelizing nodes, should be weakly positive integer
end



%% Initialize path

if option_dohpc ~= 1
    % Determine where your m-file's folder is.
    folder = fileparts(which(mfilename));
    % Add that folder plus all subfolders to the path.
    addpath(genpath(folder));
    clear folder;
end


%% Parameters

% start year
Params.t0 = 2014; % start year; used for determining time t emission intensity and for storing variables for graphing (default: 2014)

% Preference parameters
switch option_discount  % utility discount rate (annual); DICE uses 0.015 and Stern uses 0.001
    case 'nordhaus'
        Params.discount = 0.015;
    case 'stern'
        Params.discount = 0.001;
    case 'doublenordhaus'
        Params.discount = 2*0.015;
    case 'halfnordhaus'
        Params.discount = 0.5*0.015;
    otherwise
        error('Unrecognized option_discount');
end
Params.power_cons = 1.45; % coefficient of relative risk aversion (and inverse of consumption EIS) in additively separable power utility; default: 1.45; note that is not set up for log utility (eta=1)
if Params.do_ez==1
    Params.risk = 10; % coefficient of relative risk aversion; default: 10
    Params.time = 1/(1.5); % inverse of elasticity of intertemporal substitution; default: either 1/(1.5) or 1.45
end


% Consumption parameters
Params.cons_start = 55e12*1.1772; % World Bank gives year 2013 global GDP as 55 trillion in year 2005 US dollars and year 2014 US GDP deflator as 117.72
Params.cons_growth_declinerate = 0.01338; % rate at which growth rate of population declines; default: 0.01338
Params.cons_drift = @(t) exp(-3.05709)*exp(bsxfun(@times,-Params.cons_growth_declinerate,reshape((t-1+Params.t0-1960),[1 length(t)]))); % >0, a function of time since start year (make sure returns column vector of same length as t); emissions (in Gt C) per dollars of consumption (per capita in one version and absolute in another); the t-1 in the time calculation recognizes that index for 2014 is 1
Params.cons_vol = 0.013835; % default: 0.013835

% Population parameters - If set Params.pop to a constant, then population drops out entirely
Params.pop_start = 7.1e9; % World Bank gives 7.1 billion people in 2013
Params.pop_growth_declinerate = 0.01175; % rate at which growth rate of population declines; default: 0.01175
Params.pop_growth = @(t) exp(-3.82084)*exp(bsxfun(@times,-Params.pop_growth_declinerate,reshape((t-1+Params.t0-1960),[1 length(t)]))); % >0, a function of time since start year (make sure returns column vector of same length as t); emissions (in Gt C) per dollars of consumption (per capita in one version and absolute in another); the t-1 in the time calculation recognizes that index for 2014 is 1
Params.pop = @(t) Params.pop_start*exp( sum(Params.pop_growth(1:t-1)) ); % note that cannot pass a vector t in here; argument must be a scalar

% Emissions
Params.em_per_cons_declinerate = 0.014572; % rate at which emission intensity of consumption declines; default: 0.014572
Params.em_per_cons = @(t,declinerate) 0.282373*exp(bsxfun(@times,-declinerate,reshape((t-1+Params.t0-1960),[1 length(t)])))/1e12; % >0, a function of time since start year (make sure returns column vector of same length as t); emissions (in Gt C) per dollars of consumption (per capita in one version and absolute in another); the t-1 in the time calculation recognizes that index for 2014 is 1

% Damage uncertainty, from fit to Pindyck (2019)
% (note that below we set the temperature above which damages kick in)
switch option_damagetype
    case 'base'
        if option_ocean~=1
            Params.mean_logdamage_original = -6.7342; % mean of log of damage coefficient; default: -6.7342
            Params.var_logdamage_original = 1.4684^2; % variance of log of damage coefficient; default: 1.4684^2
            Params.damage_upperbound = 0.0094; % upper bound on damage parameter; default: 0.0094
        else
            Params.mean_logdamage_original = -6.5632;
            Params.var_logdamage_original = 1.4677^2;
            Params.damage_upperbound = 0.0112;
        end
    case 'level' % functional form of damages from DICE-2016R
        Params.mean_logdamage_original = -4.1653; % mean of log of damage coefficient; default: -4.1653
        Params.var_logdamage_original = 1.2826^2; % variance of log of damage coefficient; default: 1.2826^2
        Params.damage_upperbound = 0.1132; % upper bound on damage parameter; default: 0.1132
        Params.max_pctloss = 75; % maximum percentage of output that can be lost in any given year; default: 75
    case 'nordhaus' % will match expected damages to DICE-2016R
        Params.coeff_nordhaus = 0.00236; % damage coefficient from DICE-2016R
        Params.var_logdamage_original = 1.2826^2; % variance of log of damage coefficient; default: 1.2895^2        
        Params.damage_upperbound = NaN; % upper bound on damage parameter; default: NaN because irrelevant here since mean will be so small
        Params.mean_logdamage_original = log(Params.coeff_nordhaus) - 0.5*Params.var_logdamage_original; % mean of log of damage coefficient; default: log(0.00236)
        Params.max_pctloss = 75; % maximum percentage of output that can be lost in any given year; default: 75
    otherwise
        error('Unrecognized option_damagetype: %s',option_damagetype);
end
switch option_uncertain_damage
    case 'on'
        Params.var_logdamage = Params.var_logdamage_original;
        Params.mean_logdamage = Params.mean_logdamage_original;
    case 'off'
        switch option_damagetype
            case {'base','level'}
                Params.var_logdamage = 0;
                % reset location parameter so that expected damage coefficient
                % does not change with scale parameter
                Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
                    - exp(x) , Params.mean_logdamage_original);
            case 'nordhaus'
                Params.var_logdamage = 0;
                Params.mean_logdamage = log(Params.coeff_nordhaus);
        end
    otherwise
        error('Unrecognized option_uncertain_damage');
end

% Climate sensitivity uncertainty
Params.mean_logclimsens_original = 1.10704; % mean of log of climate sensitivity; default: 1.10704, from Gillingham et al. 2018
Params.var_logclimsens_original = (0.264)^2; % variance of log of climate sensitivity; default: 0.264^2, from Gillingham et al. 2018
switch option_uncertain_climsens
    case 'on'
        Params.mean_logclimsens = Params.mean_logclimsens_original;
        Params.var_logclimsens = Params.var_logclimsens_original;
    case 'off'
        Params.var_logclimsens = 0;
        % reset location parameter so that expected climate sensitivity
        % is same as under uncertainty
        Params.mean_logclimsens = Params.mean_logclimsens_original + 0.5*Params.var_logclimsens_original;
    otherwise
        error('Unrecognized option_uncertain_climsens');
end

% Carbon parameters - M1 is non-decaying stock, and M2 is geometrically decaying stock
Params.gtc_per_ppm = 2.16; % unit conversion
Params.co2_per_c = 44/12; % unit conversion
Params.M1_start = 706; % Gt C, default: 706
Params.M2_start = 148; % Gt C, default: 148
Params.M1_frac = 0.2; % fraction of emissions reaching reservoir 1; is psi_1; default: 0.2
Params.M2_frac = 0.393*(1-Params.M1_frac); % fraction of emissions reaching reservoir 2; is psi_2(1-psi_1); default: 0.393*(1-Params.M1_frac)
Params.co2_decay = 1 - 0.5^(1/300); % decay parameter for M2 reservoir; default: 1 - 0.5^(1/300)
Params.co2_preind = 280*Params.gtc_per_ppm; % in Gt C; pre-1750 CO2 cocnentration was 280 ppm (Blasing, 2014)

% Climate parameters
Params.forcing = 5.35; % change in forcing per pct change in CO2; 5.35 W/m2 is from Myhre et al. 1998
Params.seconds_per_year = 3.15576e7;
Params.inertia = Params.seconds_per_year/(0.8*10^9); % between 0 and infinity, with high values meaning no inertia and low values meaning full inertia; default: Params.seconds_per_year/(0.8*10^9)
Params.temp_vol = 0.0102; % volatility parameter in temperature process; std dev of temperature changes from 1881-2013 was 1.02% in GFDL dataset
Params.temp_start = 0.8794; % global surface temperature in 2013 was 0.8794 deg C above 1880-1899 average in NOAA dataset
if option_ocean
    Params.Cupper = 7.3; % upper ocean heat capacity; 7.3 W yr m^{-2} K^{-1} is from Geoffroy et al (2013)
    Params.Cdeep = 91; % deep ocean heat capacity; 91 W yr m^{-2} K^{-1} is from Geoffroy et al (2013)
    Params.oceanxfer = 0.74; % heat-exchange coefficent between the two layers; 0.74 W m^{-2} K^{-1} is from Geoffroy et al (2013)
    Params.oceantemp_start =  Params.temp_start*0.0068/0.85; % deep ocean temperature in 2013; default: keep ratio 0.0068/0.85 of initial ocean-to-surface temps from DICE-2016R
end

% One more damage parameter
Params.dam_temp = Params.temp_start; % temperature above which damage kicks in if option_damagetype="base" (default: Params.temp_start)


% Turn off volatilities if in that version
if ~strcmpi(option_otherloop,'forms_of_unc')
    if strcmpi(option_consvolatility,'off') && strcmpi(option_tempvolatility,'off')
        disp('Fixing both volatility parameters at zero.');
        Params.temp_vol = 0;
        Params.cons_vol = 0;
        if ~strcmpi(option_otherloop,'montecarlo_cs') && ~strcmpi(option_otherloop,'montecarlo_damage')
            Params.montecarlo_draws = 2;
        end
    elseif strcmpi(option_consvolatility,'off')
        disp('Fixing consumption volatility parameter at zero.');
        Params.cons_vol = 0;
    elseif strcmpi(option_tempvolatility,'off')
        disp('Fixing temperature volatility parameters at zero.');
        Params.temp_vol = 0;
    end
else
    if strcmpi(option_tempvolatility,'off')
        disp('Fixing temperature volatility parameters at zero.');
        Params.temp_vol = 0;
    end
end



% define vector of values will loop through for non-damage and non-drift looping parameter, with looping parameter for this run defined by option_otherloop
switch option_otherloop
    case 'off'
        otherloop_vector = 0; % placeholder
    case 'cons_vol' % baseline: 0.013835
        otherloop_vector = [0:0.005:0.04];
    case 'temp_vol' % baseline: 0.0102
        otherloop_vector = [0:0.005:0.03];
    case 'cs_var' % baseline: scale parameter of 0.264
        % loop through scale parameters; note that will adjust location parameter to hold expected climate sensitivity constant
        otherloop_vector = [0:0.1:0.3];
    case'damage_var' % baseline: scale parameter of 1.4678
        % loop through scale parameters; note that will adjust location parameter to hold the expected damage coefficient constant
        otherloop_vector = [0:0.25:1.5];
    case 'cs_and_damage_var' % loop through scale parameters for climate sensitivity and damage coefficient
        % cs scale parameter; note that will adjust location parameter to hold expected climate sensitivity constant
        otherloop_vector1 = [0:0.1:0.3];
        % damage scale parameter; note that will adjust location parameter to hold the expected damage coefficient constant
        otherloop_vector2 = [0:0.25:1.5];
        otherloop_vector = combvec(otherloop_vector1,otherloop_vector2)'; % first column is climate sensitivity scale parameter, and second column is scale parameter for damages
    case 'eta' % baseline: 1.45
        otherloop_vector = [0.5 0.99 1.01 2 2.5 3 3.5 4];
    case 'eta_with_stern' % baseline: 1.45
        Params.discount = 0.001;
        otherloop_vector = [0.5 2 3];
    case 'emint' % baseline: 0.015
        otherloop_vector = [0:Params.em_per_cons_declinerate:5*Params.em_per_cons_declinerate];
    case 'damage_exp' % baseline: 0.0014
        otherloop_vector = [0.0005 0.001 0.002]; % loop through expected values for damages; will reset location parameter accordingly if option_uncertain_damage = 'on'
    case 'damage_exp_with_stern' % baseline: 0.0014
        Params.discount = 0.001;
        otherloop_vector = [0.0005 0.001 0.002]; % loop through expected values for damages; will reset location parameter accordingly if option_uncertain_damage = 'on'
    case 'cons_growth' % baseline declines from 4.7% in 2014; will set to a constant
        otherloop_vector = [0 0.01 0.02 0.03 0.04];
    case 'discount' % baseline utility discount rate: 0.015
        otherloop_vector = [ 0.001 0.5*0.015 2*0.015 ];
    case 'montecarlo_cs'
        % Samples from climate sensitivity distribution so can see
        % how affects distribution of variables over time.  Turns off
        % consumption and temeprature volatility, and turns off uncertainty
        % about damages.
        otherloop_vector = 0; % placeholder
        Params.var_logdamage = 0;
        Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
            - exp(x), Params.mean_logdamage_original);
        Params.var_logclimsens = 0; % need later so skips quadrature
        Params.temp_vol = 0;
        Params.cons_vol = 0;
    case 'montecarlo_damage'
        % Samples from damage coefficient's distribution so can see
        % how affects distribution of variables over time.  Turns off
        % consumption and temeprature volatility, and turns off uncertainty
        % about climate sensitivity.
        otherloop_vector = 0; % placeholder
        Params.var_logclimsens = 0;
        Params.mean_logclimsens = Params.mean_logclimsens_original + 0.5*Params.var_logclimsens_original;
        Params.var_logdamage = 0; % need later so skips quadrature
        Params.temp_vol = 0;
        Params.cons_vol = 0;
    case 'forms_of_unc'
        % run determininstic model, then cons vol only, then cs uncert
        % only, then dam uncert only, then all uncert
        otherloop_vector(:,1) = [ 0; Params.cons_vol; 0; 0; Params.cons_vol]; % first col is cons vol
        otherloop_vector(:,2) = [ 0; 0; Params.var_logclimsens_original; 0; Params.var_logclimsens_original]; % second col is cs var
        otherloop_vector(:,3) = [ 0; 0; 0; Params.var_logdamage_original; Params.var_logdamage_original]; % third col is dam var        
        Params.montecarlo_draws_original = Params.montecarlo_draws;
    otherwise
        error('Unrecognized value for option_otherloop.');
end


switch option_nodamages
    case 'on'
        Params.var_logdamage = 0;
end

%%%%% End parameterization %%%%%





%% Report model setting

if Params.perturbation>0
    disp('Running perturbation analysis instead of using analytic decomposition of scc');
else
    if ~strcmp(option_damagetype,'base')
       error('Need run perturbation analysis with this damage option');
    end
end

if option_ocean
    disp('Using two-layer climate model with explicit ocean.')
end

disp(['Using ' num2str(Params.montecarlo_draws) ' Monte Carlo draws and ' num2str(Params.quadrature_nodes) ' quadrature nodes per uncertain dimension.']);
disp(['Time horizon is ' num2str(Params.horizon) ' years.  Calculating scc in first ' num2str(Params.years_for_scc) ' years.  Discretization timestep ' num2str(Params.timestep) ' years.']);
disp(['Parameter options: disc' num2str(Params.discount) ' uncert-damage-' option_uncertain_damage ...
    ' uncert-cs-' option_uncertain_climsens ]);
disp(['Other looping: ' option_otherloop]);


% make sure are not requiring log utility setup
if Params.power_cons==1
    error('Elasticity of intertemporal substitution is equal to 1, but code is not set up for log utility.');
end

% make sure timestep not greater than 1
if Params.timestep > 1
    disp('Warning: Timestep greater than 1.  Resetting to 1.');
    Params.timestep = 1;
end

if Params.do_ez==1
    if option_ocean
        error('Ocean and Epstein-Zin are incompatible options');
    end
    if ~strcmp(option_damagetype,'base')
       error('Level damages and Epstein-Zin are incompatible options') 
    end
    if Params.timestep~=1
       error('Epstein-Zin requires timestep of 1');
    end
    if Params.learn_ez == 1
        disp('Epstein-Zin, immed learning');
    else
        disp('Epstein-Zin, no learning');
    end
    Params.montecarlo_draws = 2; % irrelevant, so don't waste space with a bunch of draws
end


%% Directory setup

if option_ocean
    customdir = [customdir '_deepocean'];
end

if Params.do_ez==1
    if Params.learn_ez == 1
        customdir = [customdir '_ezprefs-learn'];
    else
        customdir = [customdir '_ezprefs-nolearn'];
    end
    customdir = [customdir '-coeffs' mat2str(Params.basiscoeffs) '-' Params.basistype '-nodes' mat2str(Params.nodes) '-' Params.nodetype];
end


if ~strcmp(option_damagetype,'base')
    customdir = [customdir '_damtype-' option_damagetype];
elseif strcmp(option_nodamages,'on')
    customdir = [customdir '_nodamage'];
end

if Params.perturbation>0
    customdir = [customdir '_perturb-' num2str(Params.perturbation)];
end

if ~strcmpi(option_otherloop,'off')
    customdir = [customdir '_loop-' option_otherloop];
end

customdir = [customdir '_vol-cons-' option_consvolatility '-temp-' option_tempvolatility];
if strcmpi(option_consvolatility,'on') || strcmpi(option_tempvolatility,'on')
    customdir = [customdir '_mc' num2str(Params.montecarlo_draws) ];
end

customdir = [customdir '_uncert-cs-' option_uncertain_climsens '-dam-' option_uncertain_damage];
if strcmpi(option_uncertain_climsens,'on') || strcmpi(option_uncertain_damage,'on')
    customdir = [customdir '_nodes' num2str(Params.quadrature_nodes)];
end

if ~strcmpi(option_otherloop,'eta') && ~strcmp(option_otherloop,'eta_with_stern') && Params.do_ez~=1
    customdir = [customdir '_eta' num2str(Params.power_cons)];
end

if ~strcmpi(option_otherloop,'discount')
    customdir = [customdir '_disc' num2str(Params.discount)];
end

customdir = [customdir '_T' num2str(Params.horizon) '_tstep' num2str(Params.timestep)];

dirsave = [dirmain '\results\' customdir];

if option_dohpc ~= 1 % create the directory only if running on personal computer
    
    if ~exist(dirsave)
        mkdir(dirsave);
    end
    disp(['Saving to ' dirsave]);
    
    cd(dirsave);
    
end


%% Loop through permutations told to run

if ~strcmpi(option_otherloop,'off')
    otherloop_steps = length(otherloop_vector);
    disp(['Doing ' num2str(otherloop_steps) ' values of parameter ' option_otherloop]);
else
    otherloop_steps = 1;
end

for index_otherloop = 1:length(otherloop_vector)
    
    tic;
    
    disp(['Beginning loop ' num2str(index_otherloop) ' out of ' num2str(length(otherloop_vector))]);
    
    % Initialize looping variable
    switch option_otherloop
        case 'off'
            % placeholder
        case 'cons_vol'
            Params.cons_vol = otherloop_vector(index_otherloop);
        case 'temp_vol'
            Params.temp_vol = otherloop_vector(index_otherloop);
        case 'cs_var'
            Params.var_logclimsens = otherloop_vector(index_otherloop).^2;
            
            % reset location parameter so that expected climate sensitivity
            % does not change with scale parameter
            Params.mean_logclimsens = Params.mean_logclimsens_original + 0.5*Params.var_logclimsens_original - 0.5*Params.var_logclimsens;
            
        case 'damage_var'
            Params.var_logdamage = otherloop_vector(index_otherloop).^2;
            
            % reset location parameter so that expected damage coefficient
            % does not change with scale parameter
            if Params.var_logdamage > 0
                Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
                    - log_normal_truncated_ab_mean( x, sqrt(Params.var_logdamage), 0, Params.damage_upperbound ), Params.mean_logdamage_original);
            else
                Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
                    - exp(x) , Params.mean_logdamage_original);
            end
            
            
        case 'cs_and_damage_var' % loop through scale parameters for climate sensitivity and damage coefficient
            
            Params.var_logclimsens = otherloop_vector(index_otherloop,1).^2;
            
            Params.var_logdamage = otherloop_vector(index_otherloop,2).^2;
            
            % reset location parameter so that expected climate sensitivity
            % does not change with scale parameter
            Params.mean_logclimsens = Params.mean_logclimsens_original + 0.5*Params.var_logclimsens_original - 0.5*Params.var_logclimsens;
            
            % reset location parameter so that expected damage coefficient
            % does not change with scale parameter
            if Params.var_logdamage > 0
                Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
                    - log_normal_truncated_ab_mean( x, sqrt(Params.var_logdamage), 0, Params.damage_upperbound ), Params.mean_logdamage_original);
            else
                Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
                    - exp(x) , Params.mean_logdamage_original);
            end
            
        case 'eta'
            Params.power_cons = otherloop_vector(index_otherloop);
        case 'emint'
            Params.em_per_cons_declinerate = otherloop_vector(index_otherloop);
        case 'eta_with_stern'
            Params.power_cons = otherloop_vector(index_otherloop);
        case {'damage_exp','damage_exp_with_stern'}
            if Params.var_logdamage > 0
                Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( x, sqrt(Params.var_logdamage), 0, Params.damage_upperbound ) ...
                    - otherloop_vector(index_otherloop), Params.mean_logdamage_original);
            else
                Params.mean_logdamage = fzero(@(x) exp(x) - otherloop_vector(index_otherloop) , Params.mean_logdamage_original);
            end
        case 'cons_growth'
            Params.cons_drift = @(t) otherloop_vector(index_otherloop); % is a constant
        case 'discount'
            Params.discount = otherloop_vector(index_otherloop);
        case 'montecarlo_cs'
        case 'montecarlo_damage'
        case 'forms_of_unc'

            Params.cons_vol = otherloop_vector(index_otherloop,1);
            Params.var_logclimsens = otherloop_vector(index_otherloop,2);
            Params.var_logdamage = otherloop_vector(index_otherloop,3);
            
            if Params.cons_vol>0
                Params.montecarlo_draws = Params.montecarlo_draws_original;
            else
                Params.montecarlo_draws = 2;
            end
                
            if Params.var_logclimsens==0
                % reset location parameter so that expected climate sensitivity
                % is same as under uncertainty
                Params.mean_logclimsens = Params.mean_logclimsens_original + 0.5*Params.var_logclimsens_original;
            else
                Params.mean_logclimsens = Params.mean_logclimsens_original;
            end
            
            if Params.var_logdamage==0
                switch option_damagetype
                    case {'base','level'}
                        % reset location parameter so that expected damage coefficient
                        % does not change with scale parameter
                        Params.mean_logdamage = fzero(@(x) log_normal_truncated_ab_mean( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), 0, Params.damage_upperbound ) ...
                            - exp(x) , Params.mean_logdamage_original);
                    case 'nordhaus'
                        Params.mean_logdamage = log(Params.coeff_nordhaus);
                end
            else
                Params.mean_logdamage = Params.mean_logdamage_original;
            end
            
        otherwise
            error('Unrecognized value for option_otherloop.');
    end
    
    
    %% Monte Carlo draws for stochasticity
    
    % Params.montecarlo_draws < 2 can cause problems once do squeeze operation
    if Params.montecarlo_draws < 2
        disp('Params.montecarlo_draws < 2 so increasing it to 2');
        Params.montecarlo_draws = 2;
    end
    
    % set up annualized covariance matrix, ordered by consumption then temperature
    covar_wiener = zeros(2,2);
    covar_wiener(1,:) = [ 1 0 ];
    covar_wiener(2,2:end) = [ 1 ];
    covar_wiener = triu(covar_wiener,0) + triu(covar_wiener,1)';
    
    % Wiener increments
    rng('default');
    for index_time = 1:Params.horizon/Params.timestep + Params.years_for_scc/Params.timestep
        % draws x processes x time
        wiener_increments(:,:,index_time) = mvnrnd( zeros(1,size(covar_wiener,1)), Params.timestep*covar_wiener, Params.montecarlo_draws);
    end
    
    % convert to Wiener paths
    wiener_paths_temp = cumsum(wiener_increments(:,:,:),3);
    % set the first entry of each path to zero and store the paths from
    % there
    wiener_paths = cat(3,zeros(size(wiener_paths_temp,1),size(wiener_paths_temp,2),1),wiener_paths_temp);
    
    % store processes at annual timestep: draws x time
    Wiener.cons = squeeze(wiener_paths(:,1,1:1/Params.timestep:end));
    Wiener.temp = squeeze(wiener_paths(:,2,1:1/Params.timestep:end));
    % store increments: draws x time
    Wiener.cons_increment = squeeze(wiener_increments(:,1,:));
    Wiener.temp_increment = squeeze(wiener_increments(:,2,:));
    
    clear wiener_paths wiener_paths_temp wiener_increments;
    
    
    
    %% Quadrature nodes for climate sensitivity
    
    % output is Params.quadrature_nodes x 1 when there is
    % uncertainty
    if isscalar(Params.var_logclimsens) && Params.var_logclimsens==0 % no uncertainty
        nodes_climsens = exp(Params.mean_logclimsens);
        weights_climsens = 1;
    else % uncertainty
        if ~strcmpi(option_otherloop,'cs_var') && ~strcmpi(option_otherloop,'cs_and_damage_var')
            [nodes_climsens, weights_climsens] = qnwlogn(Params.quadrature_nodes,Params.mean_logclimsens,Params.var_logclimsens);
        else % options: cs_var or cs_and_damage_var
            if Params.var_logclimsens>0
                [nodes_climsens, weights_climsens] = qnwlogn(Params.quadrature_nodes,Params.mean_logclimsens,Params.var_logclimsens);
            else
                nodes_climsens(:,1) = exp(Params.mean_logclimsens)*ones(Params.quadrature_nodes,1);
                weights_climsens(1,1) = 1;
                weights_climsens(2:Params.quadrature_nodes,1) = zeros(Params.quadrature_nodes-1,1);
            end
        end
    end
    
    % Reshape the warming parameter vector so
    % that the vector is along the third dimension
    Params.warming = reshape(nodes_climsens,[1 1 size(nodes_climsens,1) 1]);
    
    % convert climate sensitivity to the warming parameter
    Params.warming = Params.forcing*log(2)./Params.warming;
    
    
    
    %% Quadrature nodes for damages
    
    % Will call scripts with hard-coded nodes from Fortran90 files
    % because Matlab package seems to return incorrect nodes
    
    % output is Params.quadrature_nodes x 1 when there is
    % uncertainty
    if isscalar(Params.var_logdamage) && Params.var_logdamage==0 % no uncertainty
        nodes_damage = exp(Params.mean_logdamage);
        weights_damage = 1;
    else % uncertainty
        switch option_damagetype
            case 'base'
                
                if ~strcmpi(option_otherloop,'damage_var') && ~strcmpi(option_otherloop,'cs_and_damage_var') && ~strcmpi(option_otherloop,'damage_exp') && ~strcmpi(option_otherloop,'damage_exp_with_stern')
                    %[nodes_damage, weights_damage] = truncated_normal_rule_modified(2,Params.quadrature_nodes,Params.mean_logdamage,sqrt(Params.var_logdamage),log(Params.damage_upperbound),'empty');
                    run sub_damagenodes_base
                elseif ~strcmpi(option_otherloop,'damage_exp') && ~strcmpi(option_otherloop,'damage_exp_with_stern') % options: damage_var or cs_and_damage_var
                    if Params.var_logdamage>0
                        %[nodes_damage(:,1), weights_damage(:,1)] = truncated_normal_rule_modified(2,Params.quadrature_nodes,Params.mean_logdamage,sqrt(Params.var_logdamage),log(Params.damage_upperbound),'empty');
                        run sub_damagenodes_damvar
                    else
                        nodes_damage(:,1) = Params.mean_logdamage*ones(Params.quadrature_nodes,1);
                        weights_damage(1,1) = 1;
                        weights_damage(2:Params.quadrature_nodes,1) = zeros(Params.quadrature_nodes-1,1);
                    end
                else % option: damage_exp or damage_exp_with_stern
                    if Params.var_logdamage > 0
                        %[nodes_damage(:,1), weights_damage(:,1)] = truncated_normal_rule_modified(2,Params.quadrature_nodes,Params.mean_logdamage,sqrt(Params.var_logdamage),log(Params.damage_upperbound),'empty');
                        run sub_damagenodes_damexp % assumes scale parameter = 1.472
                    else
                        nodes_damage(:,1) = Params.mean_logdamage*ones(Params.quadrature_nodes,1);
                        weights_damage(1,1) = 1;
                        weights_damage(2:Params.quadrature_nodes,1) = zeros(Params.quadrature_nodes-1,1);
                    end
                end
                
                nodes_damage = exp(nodes_damage);
                
            case 'level'
                
                %[nodes_damage, weights_damage] = truncated_normal_rule_modified(2,Params.quadrature_nodes,Params.mean_logdamage,sqrt(Params.var_logdamage),log(Params.damage_upperbound),'empty');
                run sub_damagenodes_base
                
                nodes_damage = exp(nodes_damage);
                
            case 'nordhaus'
                
                [nodes_damage,weights_damage] = qnwlogn(Params.quadrature_nodes,Params.mean_logdamage,Params.var_logdamage);
                
            otherwise
                error('Unrecognized value option_damagetype=%s',option_damagetype);
                
        end
        
        
    end
    
    % Reshape the damage parameter vector so
    % that the vector is along the fourth dimension
    Params.cons_damage = reshape(nodes_damage,[1 1 1 size(nodes_damage,1)]);
    
    switch option_nodamages
        case 'on'
            disp('Warning: Turning damages off');
            Params.cons_damage = 0;
            weights_damage = 1;
    end
    
    
    %% Adjust if are in case where are doing Monte Carlo in either warming parameter or damage parameter (in order to explore effects on distribution of future variables)
    
    if strcmpi(option_otherloop,'montecarlo_cs') || strcmpi(option_otherloop,'montecarlo_damage')
        
        % Zero out Wiener processes
        Wiener.cons = 0*Wiener.cons;
        Wiener.temp = 0*Wiener.temp;
        Wiener.cons_increment = 0*Wiener.cons_increment;
        Wiener.temp_increment = 0*Wiener.temp_increment;
        
        if strcmpi(option_otherloop,'montecarlo_cs')
            
            seed = 123456789;
            nodes_climsens = normrnd(Params.mean_logclimsens_original,sqrt(Params.var_logclimsens_original),Params.montecarlo_draws,1);
            nodes_climsens = exp(nodes_climsens); % adjust for having been lognormal
            
            Params.warming = reshape(nodes_climsens,[length(nodes_climsens) 1 1 1]);
            
            % convert climate sensitivity to the warming parameter
            Params.warming = Params.forcing*log(2)./Params.warming;
            
        elseif strcmpi(option_otherloop,'montecarlo_damage')
            
            seed = 123456789;
            for index_sample=1:Params.montecarlo_draws
                [nodes_damage(index_sample,1), seed] = truncated_normal_b_sample( Params.mean_logdamage_original, sqrt(Params.var_logdamage_original), log(Params.damage_upperbound), seed); % need use original scale but adjusted location parameter
            end
            nodes_damage = exp(nodes_damage);
            clear index_sample seed;
            
            Params.cons_damage = reshape(nodes_damage,[length(nodes_damage) 1 1 1]);
            
        end
        
    end
    
    
    %% Divert to Epstein-Zin code if in that mode
    
    if Params.do_ez==1
        if ~isscalar(weights_damage) || ~isscalar(weights_climsens) || Params.cons_vol~=0 % there is uncertainty            
            run sub_epsteinzin;
            continue;
        else
            % in deterministic setting, Epstein-Zin differs from base model
            % only in calibration of elasticity of intertemporal
            % substitution
            Params.power_cons = Params.time;
        end
    end
    
    
    %% Initialize
    
    max_time = Params.horizon + Params.years_for_scc;
    counter_fine = 0;
    Paths.pop = zeros(1,max_time);
    if ~strcmpi(option_otherloop,'montecarlo_cs') && ~strcmpi(option_otherloop,'montecarlo_damage') % standard case
        
        Paths.M1 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths.M2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths.temp = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths.cons = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths.cons_pc = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths.forcing = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths.semielast_cons = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        
        if ~strcmpi(option_damagetype,'base')
            Paths.cons_net = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        end
        
        if option_ocean
            Paths.oceantemp = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        end
        
    else % are doing Monte Carlo draws from cs or damage distributions for illustrative purposes
        
        Paths.M1 = zeros(Params.montecarlo_draws,max_time);
        Paths.M2 = zeros(Params.montecarlo_draws,max_time);
        Paths.temp = zeros(Params.montecarlo_draws,max_time);
        Paths.cons = zeros(Params.montecarlo_draws,max_time);
        Paths.cons_pc = zeros(Params.montecarlo_draws,max_time);
        Paths.forcing = zeros(Params.montecarlo_draws,max_time);
        Paths.semielast_cons = zeros(Params.montecarlo_draws,max_time);
        
        if ~strcmpi(option_damagetype,'base')
            Paths.cons_net = zeros(Params.montecarlo_draws,max_time);
        end
    
        if option_ocean
            Paths.oceantemp = zeros(Params.montecarlo_draws,max_time);
        end
        
    end
    
    
    
    
    %% Calculate time paths
    
    for index_time = 1:max_time
        
        % Starting time
        if index_time==1
            Paths.pop(1,index_time) = Params.pop(index_time);
            if ~strcmpi(option_otherloop,'montecarlo_cs') && ~strcmpi(option_otherloop,'montecarlo_damage') % standard case
                
                Paths.M1(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.M1_start;
                Paths.M2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.M2_start;
                Paths.temp(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.temp_start;
                Paths.cons(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.cons_start;
                Paths.cons_pc(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.cons_start/Paths.pop(index_time);
                
                if ~strcmpi(option_damagetype,'base')
                    Paths.cons_net(:,index_time,:,:) = bsxfun(@times,ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1)),Params.cons_start*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,Paths.temp(:,index_time,:,:).^2))));
                    Paths.cons_pc(:,index_time,:,:) = Paths.cons_net(:,index_time,:,:)/Paths.pop(index_time);
                end
                
                if option_ocean
                    Paths.oceantemp(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.oceantemp_start;
                end
                
            else % are doing Monte Carlo draws from cs or damage distributions for illustrative purposes
                
                Paths.M1(:,index_time) = ones(Params.montecarlo_draws,1)*Params.M1_start;
                Paths.M2(:,index_time) = ones(Params.montecarlo_draws,1)*Params.M2_start;
                Paths.temp(:,index_time) = ones(Params.montecarlo_draws,1)*Params.temp_start;
                Paths.cons(:,index_time) = ones(Params.montecarlo_draws,1)*Params.cons_start;
                Paths.cons_pc(:,index_time) = ones(Params.montecarlo_draws,1)*Params.cons_start/Paths.pop(index_time);
                
                if ~strcmpi(option_damagetype,'base')
                    Paths.cons_net(:,index_time) = bsxfun(@times,ones(Params.montecarlo_draws,1),Params.cons_start*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,Paths.temp(:,index_time).^2))));
                    Paths.cons_pc(:,index_time) = Paths.cons_net(:,index_time)/Paths.pop(index_time);
                end
                
                if option_ocean
                    Paths.oceantemp(:,index_time) = ones(Params.montecarlo_draws,1)*Params.oceantemp_start;
                end
                
            end
        end
        
        % simulate stochastic processes
        if index_time < max_time
            
            M1_value = Paths.M1(:,index_time,:,:);
            M2_value = Paths.M2(:,index_time,:,:);
            temp_value = Paths.temp(:,index_time,:,:);
            cons_value = Paths.cons(:,index_time,:,:);
            if option_ocean
                oceantemp_value = Paths.oceantemp(:,index_time,:,:);
            end
            
            for index_interval = 1:1/Params.timestep
                
                counter_fine = counter_fine + 1;
                
                % Calculate changes in variables
                run sub_transitions;
                
                % Update values of variables
                M1_value = M1_value + M1change;
                M2_value = M2_value + M2change;
                temp_value = temp_value + tempchange;
                cons_value = cons_value + conschange;
                if option_ocean
                    oceantemp_value = oceantemp_value + oceantempchange;
                end
                
                
            end % end stochastic process simulation
            
            % store next period's value
            Paths.pop(1,index_time+1) = Params.pop(index_time+1);
            Paths.M1(:,index_time+1,:,:) = M1_value;
            Paths.M2(:,index_time+1,:,:) = M2_value;
            Paths.temp(:,index_time+1,:,:) = temp_value;
            if option_ocean
                Paths.oceantemp(:,index_time+1,:,:) = oceantemp_value;
            end
            switch option_damagetype
                case 'base'
                    Paths.cons(:,index_time+1,:,:) = cons_value;
                    Paths.cons_pc(:,index_time+1,:,:) = cons_value/Paths.pop(index_time+1);
                otherwise
                    Paths.cons(:,index_time+1,:,:) = cons_value;
                    Paths.cons_net(:,index_time+1,:,:) = cons_value.*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,temp_value.^2)));
                    Paths.cons_pc(:,index_time+1,:,:) = Paths.cons_net(:,index_time+1,:,:)./Paths.pop(index_time+1);
            end
            
            clear M1_value M2_value temp_value cons_value;
            clear M1change M2change forcing tempchange conschange;
            
        end
        
        % store forcing in case want it
        Paths.forcing(:,index_time,:,:) = Params.forcing*log((Paths.M1(:,index_time,:,:)+Paths.M2(:,index_time,:,:))/Params.co2_preind);
        
        % Calculations for semi-elasticity
        if ~strcmpi(option_otherloop,'montecarlo_cs') && ~strcmpi(option_otherloop,'montecarlo_damage') % standard case
            dTi_de0 = zeros(Params.montecarlo_draws,index_time,size(nodes_climsens,1),size(nodes_damage,1));
        else % are doing Monte Carlo draws from cs or damage distributions for illustrative purposes
            dTi_de0 = zeros(Params.montecarlo_draws,index_time);
        end
        for index_time_temp = 1:index_time
            
            % deriv of time index_time_temp temperature wrt time 0 CO2
            index_time_temp2 = [1:index_time_temp];
            % the various remaining permutations are just dealing with singleton dimension mismatches from the RHS to the LHS; nothing changes on the RHS
            if size(nodes_damage,1)>1
                if size(nodes_climsens,1)>1
                    adder_temp(:,index_time_temp2,:,:) = bsxfun(@times, exp(bsxfun(@plus, bsxfun(@times,-Params.inertia*Params.warming-0.5*(Params.temp_vol.^2),index_time_temp-index_time_temp2), bsxfun(@times, Params.temp_vol, bsxfun(@minus, Wiener.temp(:,index_time_temp), Wiener.temp(:,index_time_temp2)) ) ) ) ,...
                        bsxfun(@rdivide, (Params.M1_frac+Params.M2_frac*exp(-Params.co2_decay*(index_time_temp2-1))), (Paths.M1(:,index_time_temp2,:,:)+Paths.M2(:,index_time_temp2,:,:)) ) );
                    dTi_de0(:,index_time_temp,:,:) = Params.inertia*Params.forcing*trapz(adder_temp,2);
                else
                    adder_temp(:,index_time_temp2,1,:) = bsxfun(@times, exp(bsxfun(@plus, bsxfun(@times,-Params.inertia*Params.warming-0.5*(Params.temp_vol.^2),index_time_temp-index_time_temp2), bsxfun(@times, Params.temp_vol, bsxfun(@minus, Wiener.temp(:,index_time_temp), Wiener.temp(:,index_time_temp2)) ) ) ) ,...
                        bsxfun(@rdivide, (Params.M1_frac+Params.M2_frac*exp(-Params.co2_decay*(index_time_temp2-1))), (Paths.M1(:,index_time_temp2,:,:)+Paths.M2(:,index_time_temp2,:,:)) ) );
                    dTi_de0(:,index_time_temp,1,:) = Params.inertia*Params.forcing*trapz(adder_temp,2);
                end
            else
                if size(nodes_climsens,1)>1
                    adder_temp(:,index_time_temp2,:) = bsxfun(@times, exp(bsxfun(@plus, bsxfun(@times,-Params.inertia*Params.warming-0.5*(Params.temp_vol.^2),index_time_temp-index_time_temp2), bsxfun(@times, Params.temp_vol, bsxfun(@minus, Wiener.temp(:,index_time_temp), Wiener.temp(:,index_time_temp2)) ) ) ) ,...
                        bsxfun(@rdivide, (Params.M1_frac+Params.M2_frac*exp(-Params.co2_decay*(index_time_temp2-1))), (Paths.M1(:,index_time_temp2,:,:)+Paths.M2(:,index_time_temp2,:,:)) ) );
                    dTi_de0(:,index_time_temp,:) = Params.inertia*Params.forcing*trapz(adder_temp,2);
                else
                    adder_temp(:,index_time_temp2) = bsxfun(@times, exp(bsxfun(@plus, bsxfun(@times,-Params.inertia*Params.warming-0.5*(Params.temp_vol.^2),index_time_temp-index_time_temp2), bsxfun(@times, Params.temp_vol, bsxfun(@minus, Wiener.temp(:,index_time_temp), Wiener.temp(:,index_time_temp2)) ) ) ) ,...
                        bsxfun(@rdivide, (Params.M1_frac+Params.M2_frac*exp(-Params.co2_decay*(index_time_temp2-1))), (Paths.M1(:,index_time_temp2,:,:)+Paths.M2(:,index_time_temp2,:,:)) ) );
                    dTi_de0(:,index_time_temp) = Params.inertia*Params.forcing*trapz(adder_temp,2);
                end
            end
            
            clear adder_temp adder;
            
        end
        
        Paths.semielast_cons(:,index_time,:,:) = trapz( bsxfun(@times,Params.cons_damage,dTi_de0) , 2);
        
        clear dCt_dTi_divC dTi_de0;
        
    end % end time looping
    
    
    if Params.perturbation>0
        
        run sub_perturb_calcs;
        
    else
        
        
        %% Utility and scc calculations
               
        disp(['Beginning social cost of carbon calculations.'])
        
        margutil_cons_start = Params.cons_start.^(-Params.power_cons);
        
        margutil_cons_power = Paths.cons.^(-Params.power_cons); % simulation x time x climsens_node x damage_node
        
        utilloss_cons_power = Paths.cons.*margutil_cons_power; % simulation x time x climsens_node x damage_node
        
        % To get scc, sum probability-weighted components along damage dimension, 
        % then sum probability-weighted components along climate sensitivity dimension, 
        % then average along Monte Carlo simulations, 
        % then finally integrate over time
        
        % First taking expectations over damages
        scc_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),utilloss_cons_power.*Paths.semielast_cons), 4); % simulation x time x climsens_node x 1
        utilloss_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),utilloss_cons_power), 4); % simulation x time x climsens_node x 1
        semielast_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),Paths.semielast_cons), 4); % simulation x time x climsens_node x 1
        cons_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]), Paths.cons ), 4); % simulation x time x climsens_node x 1
        margutil_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),margutil_cons_power), 4); % simulation x time x climsens_node x 1
        cons_square_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]), Paths.cons.^2 ), 4); % simulation x time x climsens_node x 1
                
        % Next taking expectations over climate sensitivity ("cs" is
        % for climate sensitivity)
        scc_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), scc_dam), 3); % simulation x time x 1 x 1
        utilloss_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), utilloss_dam), 3); % simulation x time x 1 x 1
        semielast_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), semielast_dam), 3); % simulation x time x 1 x 1
        cons_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), cons_dam), 3); % simulation x time x 1 x 1
        margutil_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), margutil_dam), 3); % simulation x time x 1 x 1
        cons_square_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), cons_square_dam), 3); % simulation x time x 1 x 1        
        
        % Obtain expected values at each time ("mc" is for Monte Carlo)
        scc_mc = mean(scc_cs,1); % 1 x time x 1 x 1
        utilloss_mc = mean(utilloss_cs,1); % 1 x time x 1 x 1
        semielast_mc = mean(semielast_cs,1); % 1 x time x 1 x 1
        cons_mc = mean(cons_cs,1); % 1 x time x 1 x 1
        margutil_mc = mean(margutil_cs,1); % 1 x time x 1 x 1
        cons_square_mc = mean(cons_square_cs,1); % 1 x time x 1 x 1
        
        % Reshape to make dimensions be time x 1
        scc_mc = scc_mc'; % time x 1
        utilloss_mc = utilloss_mc'; % time x 1
        semielast_mc = semielast_mc'; % time x 1
        cons_mc = cons_mc'; % time x 1
        margutil_mc = margutil_mc'; % time x 1
        cons_square_mc = cons_square_mc'; % time x 1
        
        % loop through times in which want scc
        for scc_year = 1:Params.years_for_scc
            
            timevector = [scc_year:Params.horizon+(scc_year-1)]';
            scc_time = bsxfun(@rdivide, bsxfun(@times, scc_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), bsxfun(@power, (Paths.pop(1,timevector)'/Params.pop_start), Params.power_cons) ) ), ...
                margutil_cons_start );
            scc_nogrowthinsurance_time = bsxfun(@rdivide, bsxfun(@times, utilloss_mc(timevector,:).*semielast_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), bsxfun(@power, (Paths.pop(1,timevector)'/Params.pop_start), Params.power_cons) ) ), ...
                margutil_cons_start );
            scc_deterministic_time = bsxfun(@rdivide, bsxfun(@times, bsxfun(@power, cons_mc(timevector,:), 1-Params.power_cons).*semielast_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), bsxfun(@power, (Paths.pop(1,timevector)'/Params.pop_start), Params.power_cons) ) ), ....
                margutil_cons_start );
            scc_damagescale_time = bsxfun(@times, -Params.power_cons, bsxfun(@rdivide, bsxfun(@times, (cons_square_mc(timevector,:) - cons_mc(timevector,:).^2).*bsxfun(@power,cons_mc(timevector,:), -Params.power_cons-1).*semielast_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), bsxfun(@power, (Paths.pop(1,timevector)'/Params.pop_start), Params.power_cons) ) ), ...
                margutil_cons_start ));
            scc_precaut_time = bsxfun(@times, 0.5*Params.power_cons.*(1+Params.power_cons), bsxfun(@rdivide, bsxfun(@times, (cons_square_mc(timevector,:) - cons_mc(timevector,:).^2).*bsxfun(@power,cons_mc(timevector,:), -Params.power_cons-1).*semielast_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), bsxfun(@power, (Paths.pop(1,timevector)'/Params.pop_start), Params.power_cons) ) ), ...
                margutil_cons_start ));
            scc_growthinsur_time = bsxfun(@rdivide, bsxfun(@times, scc_mc(timevector,:) - utilloss_mc(timevector,:).*semielast_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), bsxfun(@power, (Paths.pop(1,timevector)'/Params.pop_start), Params.power_cons) ) ), ...
                margutil_cons_start );
            
            % Integrate over time
            scc(scc_year,index_otherloop) = trapz(scc_time,1);
            scc_nogrowthinsurance(scc_year,index_otherloop) = trapz(scc_nogrowthinsurance_time,1);
            scc_deterministic(scc_year,index_otherloop) = trapz(scc_deterministic_time,1);
            scc_damagescale(scc_year,index_otherloop) = trapz(scc_damagescale_time,1);
            scc_precautionary(scc_year,index_otherloop) = trapz(scc_precaut_time,1);
            scc_growthinsurance(scc_year,index_otherloop) = trapz(scc_growthinsur_time,1);
            
            % get each year's contribution to the initial year's scc
            if scc_year==1
                scc_by_year(:,index_otherloop) = scc_time*1e-9/Params.co2_per_c;
                scc_nogrowthinsurance_by_year(:,index_otherloop) = scc_nogrowthinsurance_time*1e-9/Params.co2_per_c;
                scc_deterministic_by_year(:,index_otherloop) = scc_deterministic_time*1e-9/Params.co2_per_c;
                scc_damagescale_by_year(:,index_otherloop) = scc_damagescale_time*1e-9/Params.co2_per_c;
                scc_precautionary_by_year(:,index_otherloop) = scc_precaut_time*1e-9/Params.co2_per_c;
                scc_growthinsurance_by_year(:,index_otherloop) = scc_growthinsur_time*1e-9/Params.co2_per_c;
                scc_netprecautionarydamagescale_by_year(:,index_otherloop) = scc_precautionary_by_year(:,index_otherloop) + scc_damagescale_by_year(:,index_otherloop);
                scc_residual_by_year(:,index_otherloop) = scc_by_year(:,index_otherloop) - scc_damagescale_by_year(:,index_otherloop) - scc_precautionary_by_year(:,index_otherloop) - scc_growthinsurance_by_year(:,index_otherloop) - scc_deterministic_by_year(:,index_otherloop);
            end
            
        end % end looping through scc_year
        
        % convert units from Gt C to tCO2
        scc(:,index_otherloop) = scc(:,index_otherloop)*1e-9/Params.co2_per_c;
        scc_nogrowthinsurance(:,index_otherloop) = scc_nogrowthinsurance(:,index_otherloop)*1e-9/Params.co2_per_c;
        scc_deterministic(:,index_otherloop) = scc_deterministic(:,index_otherloop)*1e-9/Params.co2_per_c;
        scc_damagescale(:,index_otherloop) = scc_damagescale(:,index_otherloop)*1e-9/Params.co2_per_c;
        scc_precautionary(:,index_otherloop) = scc_precautionary(:,index_otherloop)*1e-9/Params.co2_per_c;
        scc_growthinsurance(:,index_otherloop) = scc_growthinsurance(:,index_otherloop)*1e-9/Params.co2_per_c;
        
        % Calculate other quantities of interest
        scc_residual(:,index_otherloop) = scc(:,index_otherloop) ...
            - scc_damagescale(:,index_otherloop) - scc_precautionary(:,index_otherloop) - scc_growthinsurance(:,index_otherloop) - scc_deterministic(:,index_otherloop);
        scc_netprecautionarydamagescale(:,index_otherloop) = scc_precautionary(:,index_otherloop) + scc_damagescale(:,index_otherloop);
        
        
        % get discount rate in each year
        for index_time=1:Params.horizon
            if index_time>1
                if semielast_mc(index_time,1)>0
                    discountrate_total(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) + log( semielast_mc(index_time,:)*cons_mc(index_time,:) ) ...
                        - log( ...
                        semielast_mc(index_time,:)*cons_mc(index_time,:)^(1-Params.power_cons)/margutil_cons_start ...
                        + 0.5*Params.power_cons*(1+Params.power_cons)*(cons_square_mc(index_time,:)-cons_mc(index_time,:)^2)*semielast_mc(index_time,:)*cons_mc(index_time,:)^(-Params.power_cons-1)/margutil_cons_start ...
                        - Params.power_cons*(cons_square_mc(index_time,:)-cons_mc(index_time,:)^2)*semielast_mc(index_time,:)*cons_mc(index_time,:)^(-Params.power_cons-1)/margutil_cons_start ...
                        + (  scc_mc(index_time,:) - semielast_mc(index_time,:)*utilloss_mc(index_time,:) )/margutil_cons_start ...
                        ) );
                    discountrate_determ(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) ...
                        + Params.power_cons*log(cons_mc(index_time,:)/Params.cons_start)  );
                    discountrate_noprecaut(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) + log( semielast_mc(index_time,:)*cons_mc(index_time,:) ) ...
                        - log( ...
                        semielast_mc(index_time,:)*cons_mc(index_time,:)^(1-Params.power_cons)/margutil_cons_start ...
                        ... %+ 0.5*Params.power_cons*(1+Params.power_cons)*(cons_square_mc(index_time,index_otherloop)-cons_mc(index_time,index_otherloop)^2)*semielast_mc(index_time,index_otherloop)*cons_mc(index_time,index_otherloop)^(-Params.power_cons-1)/margutil_cons_start ...
                        - Params.power_cons*(cons_square_mc(index_time,:)-cons_mc(index_time,:)^2)*semielast_mc(index_time,:)*cons_mc(index_time,:)^(-Params.power_cons-1)/margutil_cons_start ...
                        + (  scc_mc(index_time,:) - semielast_mc(index_time,:)*utilloss_mc(index_time,:) )/margutil_cons_start ...
                        ) );
                    discountrate_noinsurance(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) + log( semielast_mc(index_time,:)*cons_mc(index_time,:) ) ...
                        - log( ...
                        semielast_mc(index_time,:)*cons_mc(index_time,:)^(1-Params.power_cons)/margutil_cons_start ...
                        + 0.5*Params.power_cons*(1+Params.power_cons)*(cons_square_mc(index_time,:)-cons_mc(index_time,:)^2)*semielast_mc(index_time,:)*cons_mc(index_time,:)^(-Params.power_cons-1)/margutil_cons_start ...
                        ... %- Params.power_cons*(cons_square_mc(index_time,index_otherloop)-cons_mc(index_time,index_otherloop)^2)*semielast_mc(index_time,index_otherloop)*cons_mc(index_time,index_otherloop)^(-Params.power_cons-1)/margutil_cons_start ...
                        ... %+ (  scc_mc(index_time,index_otherloop) - semielast_mc(index_time,index_otherloop)*utilloss_mc(index_time,index_otherloop) )/margutil_cons_start ...
                        ) );
                else
                    discountrate_total(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) ...
                        + Params.power_cons*log(cons_mc(index_time,:)/Params.cons_start) ...
                        );
                    discountrate_determ(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) ...
                        + Params.power_cons*log(cons_mc(index_time,:)/Params.cons_start)  );
                    discountrate_noprecaut(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) ...
                        + Params.power_cons*log(cons_mc(index_time,:)/Params.cons_start) ...
                        );
                    discountrate_noinsurance(index_time,index_otherloop) = Params.discount + (1/(index_time-1))*( -Params.power_cons*log(Paths.pop(1,index_time)/Params.pop_start) ...
                        + Params.power_cons*log(cons_mc(index_time,:)/Params.cons_start) ...
                        );
                end
            else
                discountrate_total(index_time,index_otherloop) = 0;
                discountrate_determ(index_time,index_otherloop) = 0;
                discountrate_noprecaut(index_time,index_otherloop) = 0;
                discountrate_noinsurance(index_time,index_otherloop) = 0;
            end
        end
        discountrate_precaut(:,index_otherloop) = discountrate_total(:,index_otherloop) - discountrate_noprecaut(:,index_otherloop);
        discountrate_riskadj(:,index_otherloop) = discountrate_total(:,index_otherloop) - discountrate_noinsurance(:,index_otherloop);
        clear discountrate_noprecaut discountrate_noinsurance;
        
        
        elapsedtime1(index_otherloop,1) = toc;
        disp(['Time spent on loop ' num2str(index_otherloop) ': ' num2str(elapsedtime1(index_otherloop,1)/60) ' minutes.']);
        
        
        
        clear scc_dam scc_cs scc_mc utilloss_dam utilloss_cs utilloss_mc semielast_dam semielast_cs semielast_mc cons_dam cons_cs cons_mc margutil_dam margutil_cs margutil_mc index_dam index_time index_cs margutil_cons_power utilloss_cons_power;
        clear timevector scc_time scc_nogrowthinsurance_time scc_deterministic_time damagescale_time scc_damagescale_time scc_precaut_time scc_growthinsur_time;
        if length(otherloop_vector)>1
            clear Paths Wiener;
            clear nodes_climsens weights_climsens nodes_climsens nodes_damage weights_damage;
        end
        
        if index_otherloop == length(otherloop_vector)
            save(['asset_results.mat']);
        else
            save(['asset_results' num2str(index_otherloop) '.mat']);            
        end
        if index_otherloop > 1
            delete(['asset_results' num2str(index_otherloop-1) '.mat']);
        end
            
        
    end
    
end % end otherloop


% Report
disp(['initial scc: ' mat2str(scc(1,:))]);
if Params.perturbation>0
    disp(['scc_exogCO2: ' mat2str(scc_exogCO2(1,:))]);
else
    disp(['initial scc_nogrowthinsurance: ' mat2str(scc_nogrowthinsurance(1,:))]);
    disp(['initial scc_deterministic: ' mat2str(scc_deterministic(1,:))]);
end



%% Save workspace

if Params.perturbation==0
    excel_scc_matrix = [ scc(1,:); scc_growthinsurance(1,:); scc_netprecautionarydamagescale(1,:); scc_precautionary(1,:); scc_damagescale(1,:); scc_deterministic(1,:) ];
else
    excel_scc_matrix = transpose([ scc_exogCO2(1,:); scc(1,:) ]);
    excel_scc_matrix(end+1,:) = excel_scc_matrix(end,:) - excel_scc_matrix(1,:);
    excel_scc_matrix(end+1,:) = 100*( excel_scc_matrix(end-1,:) - excel_scc_matrix(1,:) )./excel_scc_matrix(1,:);
end

save('asset_results.mat');



%%  LaTeX table: Create a panel for Table 1

if strcmp(option_otherloop,'forms_of_unc') && Params.perturbation==0
       
    column_labels = [ "\textbf{Source of Uncertainty}" "\textbf{Deterministic}" "\textbf{Precautionary}" "\textbf{Damage Scaling}" "\textbf{Growth Insurance}" "\textbf{Total}"];
    row_labels = [ "Consumption Volatility" "Warming Parameter" "Damage Parameter" "All Three Factors"];
    
    tablerows(1,1) = strcat("& \multicolumn{4}{c}{\textbf{Channel}} & \\ \cmidrule{2-5}  ", column_labels(1), " & ");
    for index_col=1:4
        tablerows(1,1+index_col) = strcat( num2str(column_labels(index_col+1)), "&" );
    end
    tablerows(1,6) = strcat( num2str(column_labels(end)), " \\ \midrule ");
    
    for index_row = 1:4
        tablerows(1+index_row,1) = strcat( row_labels(index_row), " & " );
        for index_col=1:5
            switch index_col
                case 1
                    tablerows(1+index_row,1+index_col) = strcat( num2str(scc_deterministic(1,index_row+1)), " & " );
                case 2
                    tablerows(1+index_row,1+index_col) = strcat( num2str(scc_precautionary(1,index_row+1)), " & " );
                case 3
                    tablerows(1+index_row,1+index_col) = strcat( num2str(scc_damagescale(1,index_row+1)), " & " );
                case 4
                    tablerows(1+index_row,1+index_col) = strcat( num2str(scc_growthinsurance(1,index_row+1)), " & " );
                case 5
                    tablerows(1+index_row,1+index_col) = strcat( num2str(scc(1,index_row+1)), " \\ " );
                    if index_row==3
                        tablerows(1+index_row,1+index_col) = strcat( num2str(tablerows(1+index_row,1+index_col)), " \addlinespace " );
                    end
            end
        end
    end
    
    % Combine into final table
    for index1=1:size(tablerows,1)
        tablefinal(index1,1) = strjoin(tablerows(index1,:));
    end
    tablefinal = strjoin(tablefinal);    
    
    save('asset_results.mat');
   
end
